%% Main script for empirical application of Goldin and Reck, "Optimal Defaults and Normative Ambiguity"

% This script generates Figures 2-4 in the final manuscript, and all
% Appendix Figures
% Our code uses structural estimates and code from Bernheim, Fradkin and Popov
% https://www.aeaweb.org/articles?id=10.1257/aer.20130907


clear all
clc

%cd to where your code is located
%cd  'XXX/finalcode/'


%Our (Goldin and Reck's) main addition is the pi parameter:
%Values of pi.
%FOR FIGURES 2a and 2b, EV BY D FOR 5 VALUES OF PI
%pivalues=(0:0.2:1);
%FOR FIGURE 2c, EV OF ACTIVE AND 6 FOR MANY VALUES OF PI
pivalues=(0:0.01:1);


matchon = 1; %Toggles the employer match on (1) and off (0).
%note: Match off no longer used.
% Company specific parameters
ms=[1,.5,.5]*matchon;

%FOR ANYTHING BUT FIGURE 2b
maxrates=[.15,.15,.15];
maxmatch=0.06;
%FOR FIGURE 2b, see what happens with a very high maximum contribution:
%maxrates=[.8, .8, .8];



%SELECT COMPANY
%We use company 3 in the paper.
%Toggle this to 1 or 2 for the results for other companies in the Appendix

company=3;

%% Parameters from BFP (Basic Model, No Anchor)

alpha=0.133969807;
t=0.2; %t in the EV and utility_optout/utility_optin functions refers to tax rate
lambda=11.81184965;
lambda2=0.401086613;
mu1=0.214986168;
mu2=0.131336922;
mu3=0.157027598;
sd1=0.090980052;

mus=[mu1,mu2,mu3]; 

%% EV Calculations

%% CALCULATIONS FOR FIG 2

% Loop over each company


%can loop over companies
%run only for company 3 for main run.
for co=company:company
% Initialize remaining parameters    
m=ms(co);
mu=mus(co);
maxcont=maxrates(co);
lowval=0.01;

%cost distribution
cost=lambda;
topcost=icdf('exp',.9999,1/cost); % Maximum value of cost for numerical integration
topnorm=icdf('norm',.9999,mu,sd1);
%topnorm = 0.99; % Maximum value of rho for numerical integration


defvalues=(0:.01:maxcont); 

dstar=zeros(length(pivalues),1);

EV=zeros(length(pivalues),length(defvalues)+1);

count=1;
% Loop over each pi value, default
for i=1:length(pivalues)
    pi=pivalues(i);
 
    for j=1:length(defvalues)
    d=defvalues(j);

    % EV
    %switchers
    temp_ev=dblquad(@(gamma,rho)ev(gamma,alpha,rho,t,m,maxmatch,d,maxcont,pi).*(1-lambda2).*cost.*exp(-cost*gamma).*normpdf(rho,mu,sd1),0,topcost,lowval,topnorm);
    temp_ev=temp_ev+normcdf(lowval,mu,sd1)*quad(@(gamma)ev(gamma,alpha,lowval,t,m,maxmatch,d,maxcont,pi).*(1-lambda2).*cost.*exp(-cost*gamma),0,topcost);
   
    %add no-cost people
    temp_ev=temp_ev+quad(@(rho)ev(0,alpha,rho,t,m,maxmatch,d,maxcont,pi).*(lambda2).*normpdf(rho,mu,sd1),lowval,topnorm);
    EV(i,j)=temp_ev;
    
    count=count+1;
        
    end
    %EV of active choice policy
    temp_ev_a=dblquad(@(gamma,rho)ev_active(gamma,alpha,rho,t,m,maxmatch,d,maxcont,pi).*(1-lambda2).*cost.*exp(-cost*gamma).*normpdf(rho,mu,sd1),0,topcost,lowval,topnorm);
    temp_ev_a=temp_ev_a+normcdf(lowval,mu,sd1)*quad(@(gamma)ev_active(gamma,alpha,lowval,t,m,maxmatch,d,maxcont,pi).*(1-lambda2).*cost.*exp(-cost*gamma),0,topcost);
   
    %add no cost people
    temp_ev_a=temp_ev_a+quad(@(rho)ev_active(0,alpha,rho,t,m,maxmatch,d,maxcont,pi).*(lambda2).*normpdf(rho,mu,sd1),lowval,topnorm);
    
    EV(i,length(defvalues)+1)=temp_ev_a;
    
    %optimal policy for given pi
    [~,istar]=max(EV(i,:));
    dtemp=[defvalues, 9.99];
    dstar(i)=dtemp(istar);
    
    %optimal policy without active choice
    [~,istar2]=max(EV(i,1:length(defvalues)));
    dstar2(i)=defvalues(istar2);
    
end

%% Results
t1=[co, pivalues]';
t2=[0.06 999 co; EV(:,7), EV(:,length(defvalues)+1), dstar];
%summarize results
results=[t1, t2]

%Plot EV over default values, not including active choice, for each value
%of pi. Toggle these on to generate the paper's figures (making sure to set
%pi to the proper set of values).

%% BUILD FIGURES 2:
%FIGURE 2a OR 2b (depending on maxcont, see above)
%{
figure('Color',[1 1 1]);
hold on
set(gca, 'fontsize',11);
plot(defvalues, EV(:,1:length(defvalues)),'LineWidth',2)
legend('pi=0','pi=0.2','pi=0.4','pi=0.6','pi=0.8','pi=1.0',...
    'Location','southoutside','Orientation','horozontal');
%Fig 2a
xticks(defvalues);
%Fig 2b
%xticks(0:0.1:0.8);
xlabel('d');
ylabel('Welfare(EV)');
%Fig 2a
ylim([-0.02, 0])
%Fig 2b
%ylim([-0.03, 0])
savefig('output/fig2b_evbydefault_co3')
print('output/fig2b_evbydefault_co3','-djpeg')
hold off
%}

%FIGURE 2c
%Compare 6% and Active choice

figure('Color',[1 1 1]);
hold on
set(gca, 'fontsize',14);
plot(pivalues, EV(:,7), pivalues, EV(:,length(defvalues)+1), 'LineWidth',2)
legend('d=0.06','d=active','Location','southoutside','Orientation','horizontal');
xlabel('pi');
xticks(0:0.1:1)
ylabel('Welfare(EV)');
yticks(-0.045:0.005:0)
savefig('output/fig2c_dstarbypi_co3')
print('output/fig2c_dstarbypi_co3','-djpeg')
%}

end
%}


%% FIGURE 4: INTERNALITIES
%Note: this takes a few minutes to run.

%values for marginal internality
mmvalues=0:0.001:0.6;
pivalues=[0,0.1,1];

%Fig 3: optimal default by pi and m:

%Now loop over pivalues and mmvalues

for co=company:company
    
    m=ms(co);
    mu=mus(co);
    maxcont=maxrates(co);
    lowval=0.01;

    %cost distribution
    cost=lambda;
    topcost=icdf('exp',.9999,1/cost); % Maximum value of cost for numerical integration
    topnorm=icdf('norm',.9999,mu,sd1);
    %topnorm = 0.99; % Maximum value of rho for numerical integration

    defvalues=(0:.01:maxcont); 

    dstar=zeros(length(pivalues),length(mmvalues));
    dstar2=zeros(length(pivalues),length(mmvalues));
    
    EV=zeros(length(pivalues),length(defvalues)+1,length(mmvalues));

count=1;

for i=1:length(pivalues)
    pi=pivalues(i);
    for k=1:length(mmvalues)
        
    for j=1:length(defvalues)
    
    d=defvalues(j);
    mm=mmvalues(k);
    
    % EV  
    %switchers
    temp_ev=dblquad(@(gamma,rho)ev_int(gamma,alpha,rho,t,m,maxmatch,d,maxcont,pi,mm).*(1-lambda2).*cost.*exp(-cost*gamma).*normpdf(rho,mu,sd1),0,topcost,lowval,topnorm);
    temp_ev=temp_ev+normcdf(lowval,mu,sd1)*quad(@(gamma)ev_int(gamma,alpha,lowval,t,m,maxmatch,d,maxcont,pi,mm).*(1-lambda2).*cost.*exp(-cost*gamma),0,topcost);
   
    %add no-cost people
    temp_ev=temp_ev+quad(@(rho)ev_int(0,alpha,rho,t,m,maxmatch,d,maxcont,pi,mm).*(lambda2).*normpdf(rho,mu,sd1),lowval,topnorm);
    EV(i,j,k)=temp_ev;
    
    count=count+1;
        
    end
    %EV of active choice policy
    temp_ev_a=dblquad(@(gamma,rho)ev_active_int(gamma,alpha,rho,t,m,maxmatch,d,maxcont,pi,mm).*(1-lambda2).*cost.*exp(-cost*gamma).*normpdf(rho,mu,sd1),0,topcost,lowval,topnorm);
    temp_ev_a=temp_ev_a+normcdf(lowval,mu,sd1)*quad(@(gamma)ev_active_int(gamma,alpha,lowval,t,m,maxmatch,d,maxcont,pi,mm).*(1-lambda2).*cost.*exp(-cost*gamma),0,topcost);
   
    %add no cost people
    temp_ev_a=temp_ev_a+quad(@(rho)ev_active_int(0,alpha,rho,t,m,maxmatch,d,maxcont,pi,mm).*(lambda2).*normpdf(rho,mu,sd1),lowval,topnorm);
    
    EV(i,length(defvalues)+1,k)=temp_ev_a;
    
    %optimal policy for given pi
    [~,istar]=max(EV(i,:,k));
    dtemp=[defvalues, 0.165];
    dstar(i,k)=dtemp(istar);
    
    %optimal policy without active choice
    [~,istar2]=max(EV(i,1:length(defvalues),k));
    dstar2(i,k)=defvalues(istar2);
    end
    
end


end

%results
dstar;


save('output/internality1.mat')

%for visualization:
dstar(1,:)=dstar(1,:)+0.001;
dstar(2,:)=dstar(2,:)+0.0005;
str='active';
    
%Take out a few instances where a computational integration error caused
%an error in the optimal default
badcols=[163,168,172,173, 174,200, 308, 382, 383, 384, 385]; 
dstar(:,badcols)=[];
mmvalues(badcols)=[];


%FIGURE 3a:
figure('Color',[1 1 1],'units','normalized','position',[.1 .1 .4 .4]);
hold on
set(gca, 'fontsize',14);
plot(mmvalues, dstar)
legend('pi=0','pi=0.1','pi=1',...
    'Location','southoutside','Orientation','horozontal');
xticks(0:0.1:max(mmvalues));
yticks(defvalues);
xlabel('marginal internality');
ylabel('optimal default');
ylim([0 0.17]);
set(findall(gca, 'Type', 'Line'),'LineWidth',2);
axes('Position',[0 0 1 1],'Visible','off');
text(0.09, 0.92, str,'fontsize',14);
savefig('output/fig3a_dstarbymm_co3')
print('output/fig3a_dstarbymm_co3','-djpeg')
hold off


%}


%Fig 3b: calculate threshold for m at which active choice no longer wins
pivalues=0:0.01:0.1;
mmvalues=0:0.001:0.3;


for co=company:company
    
    m=ms(co);
    mu=mus(co);
    maxcont=maxrates(co);
    lowval=0.01;

    %cost distribution
    cost=lambda;
    topcost=icdf('exp',.9999,1/cost); % Maximum value of cost for numerical integration
    topnorm=icdf('norm',.9999,mu,sd1);
    %topnorm = 0.99; % Maximum value of rho for numerical integration

    defvalues=(0:.01:maxcont); 

    dstar=zeros(length(pivalues),length(mmvalues));
    dstar2=zeros(length(pivalues),length(mmvalues));
    
    EV=zeros(length(pivalues),length(defvalues)+1,length(mmvalues));

count=1;
for i=1:length(pivalues)
    pi=pivalues(i);
    for k=1:length(mmvalues)
        
    for j=1:length(defvalues)
    
    d=defvalues(j);
    mm=mmvalues(k);
    
    % EV  
    %switchers
    temp_ev=dblquad(@(gamma,rho)ev_int(gamma,alpha,rho,t,m,maxmatch,d,maxcont,pi,mm).*(1-lambda2).*cost.*exp(-cost*gamma).*normpdf(rho,mu,sd1),0,topcost,lowval,topnorm);
    temp_ev=temp_ev+normcdf(lowval,mu,sd1)*quad(@(gamma)ev_int(gamma,alpha,lowval,t,m,maxmatch,d,maxcont,pi,mm).*(1-lambda2).*cost.*exp(-cost*gamma),0,topcost);
   
    %add no-cost people
    temp_ev=temp_ev+quad(@(rho)ev_int(0,alpha,rho,t,m,maxmatch,d,maxcont,pi,mm).*(lambda2).*normpdf(rho,mu,sd1),lowval,topnorm);
    EV(i,j,k)=temp_ev;
    
    count=count+1;
        
    end
    %EV of active choice policy
    temp_ev_a=dblquad(@(gamma,rho)ev_active_int(gamma,alpha,rho,t,m,maxmatch,d,maxcont,pi,mm).*(1-lambda2).*cost.*exp(-cost*gamma).*normpdf(rho,mu,sd1),0,topcost,lowval,topnorm);
    temp_ev_a=temp_ev_a+normcdf(lowval,mu,sd1)*quad(@(gamma)ev_active_int(gamma,alpha,lowval,t,m,maxmatch,d,maxcont,pi,mm).*(1-lambda2).*cost.*exp(-cost*gamma),0,topcost);
   
    %add no cost people
    temp_ev_a=temp_ev_a+quad(@(rho)ev_active_int(0,alpha,rho,t,m,maxmatch,d,maxcont,pi,mm).*(lambda2).*normpdf(rho,mu,sd1),lowval,topnorm);
    
    EV(i,length(defvalues)+1,k)=temp_ev_a;
    
    %optimal policy for given pi
    [~,istar]=max(EV(i,:,k));
    dtemp=[defvalues, 0.165];
    dstar(i,k)=dtemp(istar);
    
    %optimal policy without active choice
    [~,istar2]=max(EV(i,1:length(defvalues),k));
    dstar2(i,k)=defvalues(istar2);
    end
    
end


end

threhold=zeros(length(pivalues));
md=[mmvalues;dstar];
cut=(md(2:(length(pivalues)+1),:)>=.16).*mmvalues;
maxmm=max(cut,[],2);

maxmm(length(maxmm)+1,:)=0;
pivalues(length(pivalues)+1)=max(pivalues)+0.02;


%FIGURE 4b

figure('Color',[1 1 1],'units','normalized','position',[.1 .1 .4 .4]);
hold on
set(gca, 'fontsize',14);
plot(pivalues, maxmm)
xticks(0:0.01:max(pivalues));
xlabel('pi');
yticks(0:0.005:max(mmvalues));
ylabel('maximum marginal internality for active choice');
set(findall(gca, 'Type', 'Line'),'LineWidth',2);
savefig('output/fig3b_mcutoffsactive_co3')
print('output/fig3b_mcutoffsactive_co3','-djpeg')
hold off



save('output/internality2.mat')
%}




%% Appendix FIG 4: PLAUSIBLE INTERNALITY RANGE CHECK

%Calculate distribution of x under penalty default

%% Simulate Model
rng(1234);

N=1000; 
%mintvals=0:0.02:0.8;
mintvals=[0, 0.045];

tic
values=[];

company=3;




meanx_match=zeros(length(mintvals),1);
meanx_nomatch=zeros(length(mintvals),1);

% Specify company
i=company;
m=ms(i);


mu=mus(1,i);
maxcont=maxrates(i);


% Draw a vector of rhos
rhos=normrnd(mu,sd1,N,1);

% Initialize vector of simulated contribution rates
simulated_choices=zeros(N,1);

meanx_match=[];
meanx_nomatch=[];


for k=1:length(mintvals)
values=[];
%For first run set marginal internality=0
mint=mintvals(k);

% Loop over each draw and calculate optimal savings rate
for j=1:N
    
simulated_choices(j)=xstar_int(alpha,rhos(j),t,m,maxmatch,maxcont,mint);
end

x=tabulate(simulated_choices);
x2=zeros(19,1);

for j=0:.01:.18
    k=find(x(:,1)==j);
    if(~isempty(k))
        x2(int8((j+.01)*100))=x(k,3);
    else
        if((i<6&&j<.16)||(i>5))
        x2(int8((j+.01)*100))=0;
        end
    end
end
values=vertcat(values,x2);
values=.01*values;


rates=0:.01:.18;

meanx_nomatch=vertcat(meanx_nomatch, rates*values);

%Account for match
rates_match=(rates*(1+m).*(rates<=maxmatch)+(rates+maxmatch*m).*(rates>maxmatch));
meanx_match=vertcat(meanx_match, rates_match*values);

end


%Calculate distribution of x* under penalty default, no internalities
meanx_match(1,1)
meanx_nomatch(1,1)

OUTPUT=[mintvals', meanx_nomatch, meanx_match]

toc
save('output/internality3.mat')

%FIGURE 5 (in the Appendix):
load('output/internality3.mat')
figure('Color',[1 1 1],'units','normalized','position',[.1 .1 .4 .4]);
hold on
set(gca, 'fontsize',16);
plot(mintvals, meanx_nomatch, 'LineWidth',2, 'Color',[0 0 0])
plot(mintvals, meanx_match,'LineWidth',2,'Color',[0.5 0.5 0.5], 'LineStyle','-.')
legend('excl. match','incl. match',...
    'Location','southoutside','Orientation','horozontal');
xticks(0:0.1:max(mintvals));
xlabel('marginal internality');
ylabel('mean optimal savings rate');
ylim([0.05, 0.19]);
line([0, max(mintvals)], [0.15,0.15],'Color', [0 0 0],'LineStyle','--','LineWidth',1,'HandleVisibility','off');
line([0, max(mintvals)], [0.18,0.18],'Color', [0 0 0],'LineStyle','--','LineWidth',1,'HandleVisibility','off');
%set(findall(gca, 'Type', 'Line'),'LineWidth',2);
yticks(0:0.01:0.19);
axes('Position',[0 0 1 1],'Visible','off');
text(0.845, 0.75, 'maximum','fontsize',16);
text(0.79, 0.9, 'maximum + match','fontsize',16);
%text(0.09, 0.92, str,'fontsize',14);
savefig('output/fig4_meanxstarbymm_co3')
print('output/fig4_meanxstarbymm_co3','-djpeg')
hold off



